Charge distribution and screening in layered graphene systems. 



O 

o 
o 

Q 

in 



Instituto de Ciencia 



F. Guinea 
Materiales de Madrid. CSIC. 



Cantohlanco. 28049 Madrid. Spain 



The charge distribution induced by external fields in finite stacks of graphene planes, or in semi- 
infinite graphite is considered. The interlayer electronic hybridization is described by a nearest 
neighbor hopping term, and the charge induced by the self consistent electrostatic potential is cal- 
culated within linear response (RPA). The screening properties are determined by contributions 
from inter- and intraband electronic transitions. In neutral systems, only interband transitions con- 
tribute to the charge polarizability, leading to insulating-like screening properties, and to oscillations 
in the induced charge, with a period equal to the interlayer spacing. In doped systems, we find a 
screening length equivalent to 2-3 graphene layers, superimposed to significant charge oscillations. 

PACS numbers: 73.21.-b; 73.22.-f; 73. 21. Ac; 73.20.-r; 73.20.Hb; 73.23.-b; 73.43.-f 
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I. INTRODUCTION. 

Electronically charged systems made up of a few 
graphene layers are being intensively studied^i^i^i^i^i^iii^, 
as well as the effects of charge accumulation on the sur- 
face of bulk graphit o^'^" . Hence, the study of the charge 
distribution in finite stacks of graphene layers, or in semi- 
infinite graphite, is a problem of current interest. 

It is known that screening in a single graphene layer 
shows anomalous properties, due to the vanishing of 
the density of states at the Fermi level in neutral 
grapheneiiii^. A stack of graphene planes where elec- 
trons are confined to each layer also shows unusual 
screening propertie^^S,. In addition, defects and edges can 
lead to self doping effects in a single graphene layer-li. 

The screening properties of bulk graphite were initially 
investigated describing the system as a stack of two di- 
mensional electron fluids electrostaticsUy coupled^^. The 
electronic hybridization between layers was not included 
in the model. The simplest extension that takes inter- 
layer hybridization into account includes a tight bind- 
ing hopping element between tt orbitals at Carbon atoms 
which are nearest neighbors in adjacent layers. The intro- 
duction of this term changes substantially the electronic 
structure of bulk graphite^, and also of systems contain- 
ing few graphene layer o^^'^^ . 

The charge distribution in systems under an applied 
field must be calculated self consistently. Such a calcula- 
tion, using the interlayer hopping model described above, 
has been carried out for a graphene bilayer^^. In the fol- 
lowing, we calculate the charge distribution in graphene 
stacks of arbitrary width, and in semi infinite graphite. 
We calculate the electrostatic potential self consistently, 
and assume that the induced charge can be obtained us- 
ing linear response theory. These assumptions amount to 
the Random Phase Approximation, applied to the model 
mentioned earlier. A similar calculation for a bilayer can 
be found iniS,, and it is in reasonable agreement with 
more involved numerical calculations. 

We use, as a starting point, the calculations of the 
unpperturbed electronic structure of finite graphene 
stacks, and semiinfinite graphite, reported in™. We do 



not study the effects of other interlayer hoppings, disor- 
der, or other effects related to interactions. We also do 
not consider deviations from Bernal stacking, {1212 • • • }, 
which can alter the electronic structure at the Fermi 
leveli^. 

The following section presents the model to be stud- 
ied. We discuss then screening in semiinfinite graphite, 
and we analyze a finite graphene stack next. Section Hvl 
presents the main conclusions of our work. It also con- 
tains a discussion of the limits of validity of the results 
presented here, and their relation to previous work. 



II. THE MODEL. 
A. Electrostatic effects. 

We analyze the charge distribution at the surface layers 
of a system with many graphene layers coupled electro- 
statically to an external gate. The system is schemat- 
ically shown in Fig.[T]. A potential V is applied be- 
tween the gate and the graphene stack. An electric field, 
£ ~ V/l, where I is the distance between the gate and 
the stack. The electric potential beyond the gate is as- 
sumed to be zero, so that the voltage at the gate is —V. 
Using Gauss law we can write the total charge density 
per unit area in the stack, n as en = AttS. This charge 
is distributed among the layers, ni, 712, na, • • ■ where la- 
bel 1 stands for the outermost layer. The electric field 
in the region between layers 1 and 2 is £1-2 — 47reni/(i, 
where d is the interlayer distance. This field detrmines 
the potential in layer 2. Extrapolating this procedure to 
all layers, we find that the electrostatic potential in layer 
i satisfies: 



£4 = Ej-i -I- Ane'^dy^jik 

k=l 

or, alternatively: 

e,;_i — 2ei + e,;+i = Ane^drii 



(1) 



(2) 



We include the possible effects of a finite dielectric con- 
stant, eoj into the definition of the electrostatic charge. 



FIG. 2: (Color online. Charge susceptibility required in order 
to calculate the induced charge. 
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the directions parallel to the layers, so that the parallel 
momentum transfer in the diagram in Fig. [2] vanishes, 

q = k|i - k; = 0. 



FIG. 1: (Color online). Sketch of the system studied in the 
text. The two layers closest to the gate are charged, and lead 
to the electrostatic potential shown in the lower part. 



B. Linear response approximation. 

We analyze the system within the Hartree approxima- 
tion, and we assume that the induced charge can be well 
approximated using by linear response theory (Random 
Phase Approximation). 

The induced electron density can be written, using lin- 
ear response theory, as: 

nt = ^Xi,j<^i (3) 

j 

where Xij is a static susceptibility which describes the 
charge density induced at layer i when the potential at 
layer j is ej. The calculation of these susceptibilities is 
given by diagrams like the one shown in Fig. [2]. Their 
value is: 

k||,fcj_* k'|,fcl* k||,fci k'|,fc^ 

^^1= 1^ (4) 

where the intermediate empty and occupied states are 
labelled by their moments, ky , k|| , fc^, and /i) and j 

are layer indices. The quantities a^^" ''^^ , a^" '''"^ , " ' ^ 

k' k' 

and Uj " ' ^ in eq. ([4]) are the amplitudes of the wave- 
functions on layers i and j respectively. We will only 
consider charge distributions which are homogeneous in 



C. Electronic structure. 

We describe the electronic levels of the graphene stack 
by means of a tight binding model using the tt orbitals at 
each C atom. We consider a hopping between orbitals in 
atoms which are nearest neighbors within a given layer, 
t = 2.7eV, and a hopping between atoms in adjacent 
layers which are also nearest neighbors, t± = 0.3eV. We 
study mostly the Bernal stacking, so that the hopping 
term connects half of the orbitals in a given layer with 
half of the orbitals in the two nearest layers. At long 
wavelengths and near half filling the in plane dispersion 
is well approximated by the continuum Dirac equation, 
described in terms of the Fermi velocity, v-p = 3ta/2, 
where a — 1.4Ais the distance between Carbon atoms. 
Details of the model, and of the band structure for stacks 
with different number of layers and stacking order are 
given ii>i^. 

The system has translational invariance in the direc- 
tion parallel to the layers, so that the parallel momentum, 
k|| = {kx, ky) is conserved. If the stack is infinite, the mo- 
mentum normal to the layers, is also conserved. 

In an infinite stack with only nearest neighbor coupling 
between layers, all layers are equivalent (in general, a de- 
scription of the Bernal stacking requires two inequivalent 
layers). Then, hamiltonian for each momentum decou- 
ples in a set of 2 x 2 matrices, whose entries correspond 
to Bloch states defined in the two inequivalent sublattices 
of each layer. For a given corner of the Brillouin zone, 
the hamiltoonian can be written as: 

nj ^( 2tj_coa{k±d) vpikx +iky)\ , . 

^h^k^ - {vF{kx-lky) J 

The diagonal terms in eq. (O are determined by the inter- 
layer hopping, so that 7^22 = 0, as one of the sublattices 
is decoupled from the neighboring layers. The hamilto- 
nian in eq.([5]) reduces to the Dirac equation for t± = 0. 



D. Calculation of the charge susceptibility. 




1^ 



FIG. 3: (Color online). Sketch of the inter- and intraband 
transitions whose separate contributions to the charge sus- 
ceptibility is analyzed in the text. 



The low energy eigenenergies can be approximated as: 

t;2|kj,|2 

This approximation fails for k±d ^ 7r/2, where tha bands 
show a linear dependence on k|| . 

The allowed momenta in a finite stack with N layers 
are quantized, /c" — {mi)/[d{N + \)]^^ . In addition, the 
amplitude of the wavefunctions go to zero as sin(zA;_Lc) at 
layer i from the surface. The electronic wavefunctions in 
a finite stack are described in Appendix [XI 



The response of the electrons at both inequivalent cor- 
ners of the Brillouin Zone is the same, so that we need 
to calculate the polarizability at one K point. 

The charge susceptibility includes contributions from 
transitions between the valence and conduction bands, 
interband transitions, and trensitions within the conduc- 
tion band, intraband transitions, as schematically de- 
picted in Fig.[3j. 

We consider first the interband transitions, involving 
an occcupied state in the valence band and an empty 
state in the conduction band. We neglect the contribu- 
tion from the states with fcj_d ^ 7r/2, and we use the ap- 
proximate dispersion relation given in eq.®. Then, the 
susceptibilities in eq.Q can be written an integral over 
k||, and k±_ and k'^. The electronic states in a semiinfi- 
nite stack can also be characterized by a perpendicular 
momentum, fc^, although the corresponding wavefunc- 
tions are no longer running waves. Near the surface of a 
semiinfinite stack, the amplitudes in eq.ll]) are: 

afii''^^ =Csin(z/fcid) (7) 

where C is a normalization constant. Note that, although 
one needs in principle to define the amplitude as a two 
component spinor in each layer, the low energy states 
considered here, eq.(l6]), have vanishing amplitude on the 
sublattice connected by the interlayer hopping t_L— . 
Using eq.d?!), we finally obtain: 



^tj^ r'\ /•^^")/" ,. ^H^<t>) Mj^) sin(»0O sin(j0O cos(,/)) cos(0O f'^ dk 

"^'^ ^ ^'vlJ-n/2'^Jn/2 [cos(^) + COS (</)')] 

= -%X.log(^) (8) 



where eg is a low energy cutoff to be specified later, and we write (j) = k±d and (j)' = k'j_d. In a finite stack, these 
integrals over k± and k'j_ must be replaced by sums over the quantized momenta, see section [ill Bj . The prefactor 
^mter eq.® is defined as: 



^intcr ^ J_ ^ sin(»0) sin(j^) sin(z0O sin(j^O cos(^) cos(0O 

" 7r2 7_^/2 A/2 cos{(p) — ^-^'^ 



The values of Xij*''^'' near the boundary of the stack are plotted in Fig. [4] The limiting bulk values: 



-inter _ -inter _ 1 f^' ,A /•'"^\,.. COs[n(0 - ^Q] COs(0) COs(0O 

X™+„.™-X„ - An^J_^//'^J^^, ''^ cos(0)-cos(^') ^'"^ 



are also shown. 
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III. RESULTS. 



A. Undoped semiinfinite stack. 



In a semiinfinite stack, sufficiently far from tlie surface, so that the susceptibihties Xij have converged towards their 
bulk values, eqs.© and ^ admit the solution: 



with: 



n,; 



eoe 
noe~ 



- - 2 = ine'd J2 = log ( - 



1 + cosh(K) 
2sinh(K/2) 



arctan 



sinh 



(f 



(11) 



(12) 
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FIG. 4: (Color online). Value of the prefactor in the suscep- 
tibilities near the boundary of a stack of graphene planes, see 
eq.®. The corresponding bulk values are plotted as open 
circles on the right. 



Details of the steps involved in the derivation of eq. fll 
are given in Appendix [Bl 

In metallic systems and for k — > 1, the right hand side 
in eq.dm is equal to limfcj^_o x(^-L: k|| = 0) = D{ep), 
where -D(eF) is the density of states at the Fermi level in 
a given layer. Using this result, we obtain the Thomas- 
Fermi approximation, = ATre^dD{eF) in the limit k — *■ 
0. 

For the model for undoped graphite studied here, the 
r.h.s. of cq. (fT2)l vanishes as k ^ 0, and ordinary screen- 
ing does not take place. 

The contribution of the interband transitions to the 
charge susceptibility is not enough to give a finite charge 
compressibility at long wavelengths, leading to a behav- 
ior reminiscent of that of an insulator. Note that, at 



zero doping, the transitions between occupied and empty 
states in the limit </> — </)' ^ require a finite energy of 
order t± cos((/)), except for cf) k, (j)' k, 7r/2. These transi- 
tions have vanishing weight in the integral which gives the 
charge compressibility (see Fig.[S]). Hence, the charge 
susceptibility tends to zero at zero energy and k±d — > 0. 

Note that, on the other hand, the staggered spin and 
charge susceptibilities, k_\_d ^ 0, diverge logarithmically 
at cj ^ 0^^. This divergence is maximal for k±d — tt. 

Ea.(fT2|) admits only solutions with k > if: 



-ie^dt I 



> 6 



(13) 



The scale eo in eq.® determines the region where the 
model gives a valid description of a graphene stack. At 
zero temperature, it will be determined by interlayer hop- 
pings not considered here, disorder, and lifetime broad- 
ening due to interaction effects. At finite temperatures, 
the logarithmic divergence in eq. ([8]) is cutoff by thermal 
interband excitations, st that eg ~ T . Hence, the screen- 
ing properties of the model depend on temperature. 

At the surface, the screening by interband transitions 
is modified, as the relative strength of the different sub- 
bands, as function of /c^, is modified. The contribution 
to the real space matrix elements of the susceptibility are 
given by: 

fK/2 

d(/)sin(0z)^ sin(0j)^ cos(0) (14) 

-Tr/2 

At finite dopings the density of states of each graphene 
plane in the model studied here is: 

D{ep) = 4^ (15) 

independent of the carrier concentration. Intraband tran- 
sitions modify the charge susceptibility asfc^ — > 0. The 
charge compressibility becomes finite, leading to a bulk 
screening length : 

2 lee^dt^ 
K = T — (16) 



intra 
Xij 



8t, 



5 




J 1 I , L_ 

10 20 30 

Position 



40 50 



FIG. 5: (Color online). Transitions between the valence and 
conduction bands states for {k± — fc^)c — O.OStt and k±c = 
OAn. The apameters used are up = 1 and t± =0.1. 



FIG. 6: (Color online). Charge at the surface of a semiinfi- 
nite stack, normalized to one. Only interband transitions are 
included in the calculation. 



In the following, we use the parameters: 



e 

VF 

t±d 



— « 1 
1 

5 



(17) 



Using these values, eg. p^ gives a screening lengt due 
to bulk intraband transitions of about N 2 layers. 
The screening at the surface, however, is modified and 
reduced, as shown in eg. (fT4|) . 

Numerical results, including the full dependence on po- 
sition of Xm,m+n are shown in Fig.[B]. We have chosen 
the parameter eq such that log{t±/eo) = 6. The charge 
oscillates with every second layer. This result is consis- 
tent with the logarithmic divergence of the charge and 
spin susceptibility for k±d = tt. The induced charge is 
reduced to less than 10% of the charge in the surface 
layer in about 3 — 5 layers. 




FIG. 7: (Color online). Charge distribution in systems with 
4, 10, 50 and 200 graphene layers, normalized to one. 



B. Finite stacks. 

The eqs.® and ^ can be solved in a straightforward 
way for a system with a few layers. The discrete equiv- 
alent to the reduced susceptibilities in eq.© can be de- 
fined as: 

N N 

^inter = ^ ^ sin{kldi) sm{kldj) sm{kl di) sm{kl dj) x 

n—1 n' — l 

^ cos(fc2d)cos(fc^'d) ^ 

CfcjC^„/ [cos(fc2rf) -cos(/c2'f^)] ' 



where Ci." is a normalization constant: 

N 

Cfc-^ =^sin2(fc^d) (19) 

A finite bilayer can be charged, so that the Fermi en- 
ergy does not need to lie exactly between the valence 
and conduction bands. Then, we have to add to the sus- 
ceptibility a contribution from intraband transitions at 
k|| — > 0. The contribution of each subband is determined 
by its contribution to the total density of states, and it 
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is independent of the position of the Fermi level. Hence, 
intraband transitions lead to: 



2t±_ 



E 



sin2(fc^dj) sin^ {kidj) |cos (/c^d) 



CI, 



The charge distribution is given by the equations: 



(20) 



k=l 



N 



n. = '^o + E(^^r + x'D^^ (21) 

where rto is a constant which fixes the total density, in 
turn detrmined by the applied field between the gate and 
the stack. 

The diagonal intraband susceptibilities near the sur- 
face of a semiinfinite stack satisfy: 



-intra 
Ann ^ 



tt/2 



d(f> sin* {n(j)) cos(0) 



-7r/2 



1 



(-1)" 

4(16n2 - 1) 2n2 - 1 



(22) 



showing even-odd oscillations as function of the distance 
to the surface, n, as well as a slow convergence to the 
bulk limit, n — > oo. 

If we neglect the intraband susceptibility, we find for a 
bilayer: 



ei - £2 



27re^(ino 



1 - 27re2£i(xii - X12) 



(23) 



where n is the total carrier density, and xii ^nd xi2 
are the bilayer interband susceptibilities, defined using 
eqs.lfTSl) and By choosing the low energy cutoff, eo 
such that eo = ep = (7rupno)/^i; we recover the results 
ini^ to lowest order in no- 

Examples of the charge distribution in doped systems 
with different number of layers are shown in Fig. [7]. We 
find that: 

i) The charge distribution is rather homogeneous in 
narrow stacks, and it becomes concentrated at the sur- 
faces for stacks wider than the screening length. 

ii) There are oscillations as function of the distance to 
the surface in wide stacks, as in the case of a semiinfinite 
stack. These oscillations can lead to charge with different 
sign in neighboring layers. 



IV. CONCLUSIONS. 

We have calculated the charge distribution in stacks 
of graphene layers in an applied field. We describe the 
electronic structure by a tight binding model for the tt 
orbitals, which includes hopping between sites which are 



nearest neighbors in adjacent layers. The self consis- 
tent electrostatic potential is obtained assuming linear 
response theory, so that our approximations amount to 
the Random Phase Approximation. 

The electronic structure of these systems shows a va- 
lence and conduction bands with parabolic dispersion as 
function of the parallel momentum, for all values of the 
perpendicular momentum. These bands touch at zero en- 
ergy, which defines the chemical potential in the undoped 
case. Finite and infinite stacks are gapless. 

The charge susceptibility can be written as a sum of 
intra- and interband contributions. In a clean, neutral 
system only interband transitions contribute. In this 
case, an infinite stack of layers shows non metallic screen- 
ing, as the long wavelength charge susceptibility vanishes. 

The charge susceptibility for finite momenta perpen- 
dicular to the layers shows a logarithmic divergence at 
low energies. This divergence is maximal for a wavevec- 
tor inversely proportional to the distance between the 
layers, leading to charge oscillations. In addition, the 
logarithmic divergence at finite wavevectors implies that 
the screening properties of undoped systems can show a 
depedence on the frequency at which they are probed, or 
on temperature. 

The charge distribution induced at a surface can be 
extended over many layers, and it shows a decaying mod- 
ulation with a period equal to the distance between the 
layers. 

In doped samples, the long wavelength charge polariz- 
ability is finite, and it is dominated by intraband transi- 
tions. We find that an external field is screened within 3-5 
layers from the surface. The screening length is indepen- 
dent of carrier concentration. The interband transitions 
lead to oscillations in the induced charge, as in the un- 
doped case. As these fluctuations depend logarithmically 
on a low energy cutoff, they can also show a dependence 
on temperature or frequency. 

The screening length in doped stacks obtained here 
depends on the values of the parameters given in (jl7p . 
Their bulk values are not known with precision^S, and it 
is possible that some of them, like the effective electric 
charge or the interlayer hopping, change near an inter- 
face. Hence, the value of 3-4 layers over which the charge 
is delocalized in doped samples is only approximate. 

Our calculation does not take into account effects such 
as next nearest neighbor hoppings, disorder, or deviations 
from linear reponse. The existence of other hoppings, or 
disorder, will define a low energy scale below which the 
results will no longer be valid. 

Deviations from linear response theory depend on the 
strength of the induced electrostatic potential with re- 
spect to the parameters which define the band structure, 
the smaller of which is the interlayer hopping, t^, which 
we have assumed to be t±_ ^ 0.1 — 0.3eV. Typical differ- 
ences in electrostatic potentials between adjacent layers 
are given by — ei_i « e^dn, where n is the induced 
charge per unit area. For n ~ 10^^ — 10^^ cm"^, we ob- 
tain ei — Ci-i ^ 10~^ — 10~2eV, so that the assumption of 
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linear response is probably valid. Finally, we have only 
studied the Bernal stacking, 1212 ■ • • . Regions of rhom- 
bohedral stacking reduce the density of states near the 
Fermi leveli^, and they will decrease the screening in the 
system. 

It is interesting to consider the similarities and differ- 
ences of this work with the related calculation inAS,. In 
this reference, it was assumed that a low density two di- 
mensional electron gas existed in each layer, and that the 
electrons could not move between different layers. This 
approximation is equivalent to consider only the diago- 
nal susceptibility, x„„, in eq.®. For the 2DEG, only 
intraband transitions exist, leading to Xnn = D{ep) — 
7iv X {m* /{nh'^), where is the number of valleys. This 
assumption leads to metallic screening, with a decay of 
the electrostatic potential into the bulk, eq. (|12p . given 
by « {e^dn^m*)/{TTh^). In our case, the existence of 
a finite density of states near the Fermi level is a conse- 
quence of the finite interlayer hopping, t^, as the elec- 
tronic structure reduces to a stack of decoupled Dirac 
equations in the absence of hopping. The charge oscilla- 
tions near a surface are directly related to the quantum 
coherence between stacks induced by the interlayer hop- 
ping. Hence, they cannot be obtained in the model used 
ii*i^. 
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APPENDIX A: ELECTRONIC 
WAVEFUNCTIONS IN FINITE STACKS. 

We define the amplitude of the wavefunction with 
parallel momentum ky , perpendicular momentum /c" at 
layer M as a two component spinor: 

where ag.^^ j,„ and (3^^^ refer to the amplitudes on 

the A sublattice, whose atoms in different layers are near- 
est neighbors, and the B sublattice, whose atoms in dif- 
ferent layers are not connected. In order to satisfy the 
open boundary conditions at the surfaces of the stack, 
these amplitudes must be of the form: 

( J""" ] siniklcM) (A2) 

and: 



2t± cos(fc"c) VF{kx + iky) 



(A3) 



The low energy eigenvalues, |eg. j,„ | ^ i_L are given by: 



= cos(fclc) T V*i cos2(fc^c) + .||k|| P « ^^-^-^^ (A4) 



within this low energy approximation, eq. 
that: 



a 



k||,fc^ 



C. 



1/2 



3]) implies 



(A5) 



n in eq. pni is through the factor cos[n{(j> — (/>')]. In or- 
der to regularize the summations over n, we adda small 
decaying factor, e: 



where C^i"^ is a normalization constant. 

APPENDIX B: CALCULATION OF THE BULK 
SUSCEPTIBILITY. 



cos[rt(^ — (/)')] 



e coslTi 



[nicp - 0')] (Bl) 



We derive the right hand side of eq. ([72]) using the ex- 
pressions in eq. pO|) . The dependence on the layer index 



For K = in cq. p2)) the n dependent part of the sum 
gives: 
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E 



g-eln|gin(0o+0-0') 



sinh(e) 



e-»0 



> 2TTS{(f>0 + - 0') 



n— — oo 



cosh(e) — cos(0o + — 0') 
where, for convenience, we have also incuded a shift, 0o- Using this result, we also obtain: 



(B2) 



Tr/2 fZixll 



-7r/2 J7r/2 



6' ^ cos[7i(0o + 0-0')]- 



, cos(0) cos(0') /"^/^ cos(0) cos(0o + 0) 



n— — oo 



C0S((/)) — cos(0') y_^/2 cos(0o + 0) ^ C0S((/)) 



1 + cos(0o) 
4sin(0o/2) 



log 



1 + sin(0o/2) 



1 -sin(0o/2) 



(B3) 



The summation for k 7^ in eq. (|12p is formally equivalent 
to the replacement 0o ^ zk. Making this substitution, 
we find: 



J- (in) 



1 + cosh(K) 
2sinh(K/2) 



I 

For K ^ we find: 



arctan 



sinh (!)]-! (B4) 



180 



(B5) 
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